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vo Abstract 
(N 

We investigate a computational device that harnesses the effects of Bose-Einstein condensation 
(BEC) to accelerate the speed of finding the solution of a given optimization problem. Many 
computationally difficult problems, including NP-complete problems, can be formulated as a ground 
state search problem. In a BEC, below the critical temperature, bosonic particles have a natural 
tendency to accumulate in the ground state. Furthermore, the speed of attaining this configuration 
is enhanced as a result of final state stimulation. We propose a physical device that incorporates 

O 

these basic properties of bosons into the optimization problem, such that an optimized solution is 

in 

found by a simple cooling of the physical temperature of the device. We find that the speed of 
convergence to the ground state can be sped up by a factor of N at a given error, where N is the 
boson number per site. 
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Quantum computation promises to offer great increases in speed over current computers 
due to the principle of superposition, where information can be processed in a massively 
parallel way pQ. The quantum indistinguishability [2] of particles, another fundamental 
principle of quantum mechanics, remains relatively unexplored in the context of information 
processing. Bosonic indistinguishability is the mechanism responsible for phenomena such 
as Bose- Einstein condensation (BEC) [3J. We show that by using bosonic particles it is 
possible to speed up the computation of a given optimization problem. The method takes 
advantage of the fact that bosonic particles tend to concentrate in the minimal energy state 
at low temperatures. Since many difficult computational problems can be reformulated as 
an energy minimization problem jl], this is attractive for such computational purposes that 
a large number of bosons lie in the ground state configuration. The origin of the speedup 
is due to bosonic final state stimulation, an effect that is familiar from stimulated emission 
of photons in lasers |5j. This allows the system to move towards the ground state at an 
accelerated rate. 

We formulate the computational problem to be solved as an energy minimization problem 
of an Ising Hamiltonian [4j. For example, the NP-complete MAX-CUT problem [6], where 
the task is to group M vertices into two groups A and B such as to maximize the number 
of connections between the groups, is known to be equivalent to the Hamiltonian Hp = 
JijCTiUj, where Jjj is a real symmetric matrix that specifies the connections between 
the sites i,j, and <7j = ±1 is a spin variable. The task is then to find the minimal energy 
spin configuration {<7j}. In simulated annealing [7], very long annealing times are necessary 
to ensure that the system does not get caught in local minima. Quantum annealing j8] 
overcomes such problems due to local minima by introducing a quantum tunneling term but 
requires a slow adiabatic evolution to prevent leaks into excited states. 

The computational device we have in mind is shown in Figure [TJ Each spin <jj in Hp is 
associated with a trapping site containing iV bosonic particles. The bosons can occupy one 
of two spin states, which we label by a = ±1. Any particle that displays bosonic statis- 
tics with an internal spin state may be used, such as exciton-polaritons in semiconductor 
microcavities, which have recently observed to undergo BEC [9HTT] or neutral atoms with 
an unpaired electron in atom chips [12]. Systems that undergo BEC are natural choices 
for implementation of such a device, since similar principles to the formation of a BEC 
are required in order for the rapid cooling to the solution of the computational problem. 
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FIG. 1: Each site of the Ising Hamiltonian is encoded as a trapping site, containing N bosons. The 
bosons can occupy one of two states a = ±1, depicted as either red or blue. The interaction between 
the sites may be externally induced by measuring the average spin on each site i via the detectors, 
which produce a detector current hit). A local field on each site equal to Bi = rj^- Jijlj{t)/ \ff\ 
is applied via the feedback circuit. The system dissipates energy according to the coupling a to 
the environment. 

Exciton-polaritons possess a spin of a = ±1 which can be injected by optical pumping with 
right or left circularly polarized laser beam. The sites are externally controlled such as to 
follow the Hamiltonian 

H = ^ J ij S i S j (!) 

where Si = Y^,k=i a i * s ^ ne total spin on each site i, and is the same matrix as in Hp 
which specifies the computational problem. The ground state spin configuration of ([!]) is 
equivalent to the original Ising model Hamiltonian Hp [21\. This can be seen by noting 
that the same spectrum as Hp is obtained when the site spin is maximized \Si\ = N. The 
energy between these levels connect linearly as the spin on a particular site is changed from 
Si = —N to N or vice versa. 

The interaction Hamiltonian Q may be produced by measuring the total spin Si on 
each site, processing those measurement results and feeding an appropriate control signal 
back into the system by applying a local dc field on site i. For example, say at a particular 
instant a spin measurement of all the sites are made, giving the result {Sj}. Then at that 
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moment a local field Bi = ^) . JijSj is applied on site i, yielding the effective Hamiltonian 
H = Y^iBiSi. The measurement and the feedback process are continuous. Although Ju- 
lias a large connectivity and is long-ranged, by using such a feedback method to induce the 
interactions there is no restriction to the kind of interactions that can be produced in 
principle. The above argument can be formulated in the framework of quantum feedback 
control. We start with the Wiseman-Milburn feedback master equation [15] dp c /dt = C^p c + 
D[C]p c — iy/rj[F, Aip c ] + D[F]p c , where £ is a Liouville superoperator describing the internal 
dynamics of the system, D[C]p = CpC^ — {C^C, p}/2 is the Lindblad superoperator, C is 
the measurement operator due to the meter coupling, rj is the detector efficiency, Ai is the 
measurement superoperator, and p c is the density matrix of the system conditional on prior 
measurement outcomes. We consider Markovian feedback and the system is acted on by 
a Hamiltonian i^ fb (t) = I(t)F, where I(t) is the feedback current due to the measurement 
outcome [T6] . 

We now define each of the variables in the master equation for our specific implementation. 
Our system consists of a set of cross-coupled systems such as that shown in Fig. [TJ First 
consider one particular site i. The meter measures the z-component of the spin, thus we have 
C = \J^jS* = a/7(— 2^j_ + N), where 7 is the rate constant representing the measurement 
strength, nj_ is the number operator counting the number of down spins on site i, and 
we have assumed N bosons per site. In order that the system can dissipate energy out of 
the system we have a dissipation term on each site C Q p c = aD[S^]p c , where a is a rate 
constant determining the time scale of the dissipation (cooling), = a\__a i+ , and a ia is the 
annihilation operator for a boson on site i in the state a — ±1. The first two terms of the 
master equation thus describe a cooling process with a dephasing term originating from the 
measurement of the ^-component of the spin. The back-action of the ^-measurement gives 
a measurement superoperator M.p = S-p + pSf. As a result of the feedback, on each site 
we apply a field in the z-direction such that F oc Sf. 

Now consider the complete feedback system as a whole. Consider applying a feedback 
Hamiltonian of the form H ib (t) = T SiJijIj(t)/-\/Vi where Ij(t) is the current resulting 

from the measurement of site j, Jij is the same matrix specifying the problem Hamiltonian 
([T]), and T is a overall constant. Inserting these expressions into the feedback master equation 
gives dp/dt = J2i [aD[Sr]p + lD [St]p + g J%D[St\p - 1 Ejd rj ^ S jP + P S i{ • 
Due to the symmetric nature of the matrix, the last term in the above equation can be 
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written [19j as —iT[H,p], where H is given in equation Q. This gives the time evolution 
of the density matrix dp/dt = -iT[H, p] + a EiD[Sr]p + E,(g 4 + l)D[S*\p. The 
first term is an evolution of the system according to the Hamiltonian ([I]), which shows that 
the feedback Hamiltonian H^t) indeed reproduces the desired Hamiltonian ([!]). The second 
term is a cooling of the system as before, and the third is a dephasing term originating from 
the measurement on each site, as well as a contribution from the feedback circuit noise. 

Initially each site is prepared with equal populations of a = ±1 spins, which can be 
achieved by using a linearly polarized pump laser, in the case of exciton-polaritons. The 
system is cooled in the presence of the interactions between the sites, by immersing the 
system in an external heat bath. The readout of the computation is simply performed by 
measuring the total spin on each site after the system cools down by dissipating heat into 
the environment. The sign of the total spin gives the information of <Ji = ±1 for the original 
spin model. Since the "computation" here is the cooling process itself, no complicated gate 
sequence needs to be employed to obtain the ground state. 

To understand the effect of using bosons, first compare the thermal equilibrium config- 
uration of a system described above with an equivalent system that uses classical, dis- 
tinguishable particles. As a simple example, consider the two site Hamiltonian H = 
— JS1S2 — XN(Si + S2), where the second term is included such as there is a unique ground 
state in spite of the Si <H- — Si symmetry of the first term in the Hamiltonian. For a single 
spin on each site and J, A > 0, the ground state configuration is a% = 1, a 2 = 1, which 
we regard as the "solution" of the computational problem. We neglect the presence of an 
on-site particle interaction cx Sf here since we assume that the strength of the interactions 
J produced by the induced feedback method can be made much larger than such a term 
which may occur naturally. 

In Figure 2 we show the average spin on a single site of the two site Hamiltonian, which 
can be calculated from standard partition function methods accounting for bosonic counting 
factors. Comparing bosonic particles and classical distinguishable particles, we see that the 
bosonic case has a larger average spin for N > 1 and all temperatures, corresponding to 
a spin configuration closer to the ground state. As the particle number is increased, the 
temperature required to reach a particular (Si) increases. For the bosonic case, the required 
temperature increases linearly with N, while for distinguishable particles it behaves as a 
constant for large N. This results in an improved signal to noise ratio for the bosons in 
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FIG. 2: The average spin of the two site Ising Hamiltonian as a function of the boson number iV 
and rescaled temperature ksT/JN. (a) indistinguishable bosons and (b) classical distinguishable 
particles. The parameters used are J = 10 and A = 0.5. 

comparison to distinguishable particles. The concentration of particles in the ground state 
configuration for bosons is precisely the same effect that is responsible for the formation of 
a BEC. Since the ground state corresponds to the solution of the computational problem, 
this corresponds to an enhanced probability of obtaining the correct answer at thermal 
equilibrium. 

We now turn to the time taken to reach thermal equilibrium, after initially preparing the 
system with equal populations of a = ±1 particles on each site. We generalize the methods 
of Ref. [32] to our bosonic Ising model, also accounting for transitions beyond first order in 
perturbation theory. Given the M-site Hamiltonian H = JijSiSj + XN £\ Si, the states 
are labeled \k) = YlfLi i 1 = (a\ + ) ki (a\_) N ~ ki \0), where the range from to N, a\ a is 

v/ hi !(7V— hi) ! 

the creation operator for a boson on site % in the state a, and we have defined the vector 
k = (hi, k 2 , . . . , fcjif). The probability distribution then evolves according to 

dp M N ~ ki 

where 5ki = (0, . . . , 0, Ski, 0, . . . , 0) is a vector in the direction of the iih axis of the M- 
dimensional hypercube. The w(k,5ki) is a weight factor for the process \k) — > \k + Ski), 
containing a transition rate factor from Fermi's golden rule and a (1 ± 7) factor to ensure 
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FIG. 3: (a) Equilibration of a two level system with energy levels E\ = and Ei = 10 at ksT = 10 



and boson numbers as shown, (b) The equilibration time for the 4 site Ising model with Ji 
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and A = — 1. (c) The 1/N dependence of the 4 site Ising model in (b). (d) The residual energy 
for the 2 site Ising model with J = 10 and A = 0.5 after annealing the system with an exponential 
schedule with time constant to for various boson numbers as shown. All calculations use a = 1 
and f = 0.001. 

that the system evolves to the correct thermal equilibrium distribution, in a similar way to 
that discussed in Ref. [13]. We calculate the weight factors to have the form 



w(k,5ki) = (l + ^(5k))- 



((6k- l)!) 2 * 

where a is a rate constant determining the overall timescale, 
tanh[-<$fy3 (A AT + J^(2% - N)\ and 



(2) 



H(8k) 



]\ & * =l (k + m)(N - k - 5k + m) 5k > 



F(k,5k) 



rim=# - 1^1 +m)(N -k + m) 5k < 



The J^/c, <$fc) factors are final state stimulation factors due to bosonic statistics, from ma- 
trix elements \(k + (5A: | V^' 5 ^ | A:) | 2 in Fermi's golden rule, where the perturbation causing the 
transition is V oc a+a_ + a)_a + [20]. Transition beyond order one are suppressed by the 
coefficient £ <C 1. 

We use the standard numerical differential equation solver supplied by Mathematica to 
evolve pk for small boson numbers. Figure [3k shows the cooling of the system for A^ = 1 



and N = 50 particles. As the number of particles is increased, we see that the time taken to 
reach equilibrium is considerably reduced, as well as a high proportion of particles occupying 
the ground state. For low temperatures, the time dependence of the single site case can be 
approximated by the rate equations ^ = — = a(ni + l)n2, where are the populations 
on levels i — 1,2. Analytically solving this gives a equilibration time of r ~ 1/ocN for large 
N, explicitly showing the final state stimulation effect. 

As is well known from past studies of simulated annealing, the presence of local minima 
slows down the time for equilibration dramatically. To illustrate the behavior of the system 
in the presence of local minima, we choose a typical four site Ising model with a local 
minimum state 4-4-4—4- an d a global minimum state tttt- We define the error probability e 
as the probability of failing to obtain the correct ground state configuration after a single 
measurement of the total spin: e = 1 — ^exp[— H/ksT]/Z , where Z is the partition function 
and the summation is over all configurations with the same sign of spin as the correct ground 
state. The effect of the local minimum can be seen from the N = 1 case shown in Figure 
[3b, where the time rapidly increases as e — > (i.e. T — » 0) unlike the single site case. In our 
simulations, we assume that the Hamiltonian is correctly implemented by the feedback 
scheme, and use a kinetic Monte Carlo method [22] to numerically calculate the cooling time 
starting from a T = oo configuration. A final thermal equilibrium temperature is set, which 
determines the error probability. Fig. ^jp shows that as the boson number is increased, there 
is a significant speedup at constant error of several orders of magnitude. There is a small 
odd/even effect due to the definition of the error e. The curves approach zero equilibriation 
time as ex 1/N for large N (Fig. [3^). In all our numerical simulations we have found that 
bosons are able to speed up the equilibration times, with a oc 1/N speedup for large N, in 
a similar way to the single site example. 

The scheme is also compatible with a thermal annealing procedure, where the temperature 
is gradually reduced to zero starting from a high temperature configuration. We calculate 
the residual energy, defined as the average energy above the ground state of the system 
following the annealing procedure. An exponential annealing schedule with time constant 
To is used, starting from a temperature corresponding to an error of e = 0.7. Times up to 
4r are annealed where the system no longer responds to the cooling. Fig. [3ji shows that 
the residual energy is suppressed for all N > 1, thus again displaying an improvement due 
to bosonic final state stimulation. 
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We conclude that the scheme as shown in Figure [T] has a systematic way of improving 
the standard Ising model, in terms of a speedup proportional to the number of bosons N 
per site. The origin of the speedup can be understood in the following simple way. The use 
of many bosons increases the energy scale of the Hamiltonian from ~ to ~ NJ^. Due 
to bosonic statistics, the coupling of the spins to the environment is increased by a factor of 
~ N. Thus by constructing a system out of bosons we have increased the energy scale of the 
entire problem by a factor of N, which results in a speedup of N. Spin flips due to random 
thermal fluctuations also occur on a timescale that is faster by a factor of N, resulting in 
a faster escape time out of local minima. We emphasize that although the device discussed 
in this letter is a computational device that uses quantum effects, it is rather different to a 
quantum computer, since the off-diagonal density matrix elements of the state of the device 
are explicitly zero at all times. For these reasons we expect the scaling of the equilibration 
time with the site number M is not faster than exponential, in analogy to the classical case. 
The speedup then manifests itself as a suppressed prefactor of this exponential function, 
which can be accelerated by a factor of N. In its present form, the device can simulate any 
kind of optimization problem that can be written as an Ising model involving two spins, 
such as the graph partitioning problem, 2SAT, MAX-2SAT, and others. Extension of the 
device to involve fc-body interactions give a natural implementation of problems such as 
fc-SAT (k > 3). 
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